library(tidyverse)
library(devtools)
library(nnet)
library(lubridate)
library(eeptools)
library(data.table)
library(radiant.data)
library(sjPlot)
library(dplyr)
library(rctutils)
library(mlogit)
library(lme4)
library(ggplot2)
library(stringr)
library(MASS)
library(ggmap)
library(zipcodeR)
library(brant)
library(logistf)
library(car)


cdr <- read_csv('/Users/celinascott-buechler/DFP/cdr/cdr_cleaned.csv', col_names =TRUE)
pp <- read_csv('/Users/celinascott-buechler/DFP/cdr/Power_Plants.csv', col_names=TRUE)

#pp_full <- left_join(pp, cdr, by = c("Zip" = "ZipCode"),relationship = "many-to-many")
#pp$zip_rad <- search_radius(pp$Y, pp$X, radius = 5)

pp_petro <- pp %>% filter(PrimSource=='petroleum')
pp_ng <- pp %>% filter(PrimSource=='natural gas')
pp_coal <- pp %>% filter(PrimSource=='coal')

pp_petro$petro_bin <- 1
pp_petro <- pp_petro %>% group_by(Zip) %>% summarise(petro_sum = n(), petro_capacity_mean = mean(Total_MW), petro_capacity_sum = sum(Total_MW))
pp_ng$ng_bin <- 1
pp_ng<- pp_ng %>% group_by(Zip) %>% summarise(ng_sum = n(), ng_capacity_mean = mean(Total_MW), ng_capacity_sum = sum(Total_MW))
pp_coal$coal_bin <- 1
pp_coal<- pp_coal %>% group_by(Zip) %>% summarise(coal_sum = n(), coal_capacity_mean = mean(Total_MW), coal_capacity_sum = sum(Total_MW))

cdr <- left_join(cdr, pp_petro, by = c("ZipCode" = "Zip"), relationship = "many-to-many")
cdr$petro_sum[is.na(cdr$petro_sum)] <- 0
cdr <- left_join(cdr, pp_ng, by = c("ZipCode" = "Zip"), relationship = "many-to-many")
cdr$ng_sum[is.na(cdr$ng_sum)] <- 0
cdr <- left_join(cdr, pp_coal, by = c("ZipCode" = "Zip"), relationship = "many-to-many")
cdr$coal_sum[is.na(cdr$coal_sum)] <- 0

cdr$ff_sum <- cdr$petro_sum+cdr$ng_sum+cdr$coal_sum

cdr$petro_capacity_sum[is.na(cdr$petro_capacity_sum)] <- 0
cdr$ng_capacity_sum[is.na(cdr$ng_capacity_sum)] <- 0
cdr$coal_capacity_sum[is.na(cdr$coal_capacity_sum)] <- 0

cdr$ff_cap_sum <- cdr$petro_capacity_sum +cdr$coal_capacity_sum +cdr$ng_capacity_sum

#set reference levels ahead of running statistical tests
cdr <- relevel_columns(cdr, party3 = "Democrat", Education = "No high school diploma", Income = "Under $25,000", Race= 'White', education3="High school or less")

table(cdr$Gender)/length(cdr$Gender)*100
table(cdr$Race)/length(cdr$Race)*100
table(cdr$Education)/length(cdr$Education)*100
table(cdr$Party)/length(cdr$Party)*100

cdr$party3[cdr$party3=="Neither"] <- "Independent"


####Messenger####

cdr_messenger <- cdr %>%
  mutate(support = case_when(
    !is.na(DemMessenger) ~ DemMessenger,
    !is.na(RepMessenger) ~ RepMessenger,
    !is.na(BipartMessenger) ~ BipartMessenger,
    TRUE ~ NA_character_
  ))

cdr_messenger <- cdr_messenger %>%
  mutate(support2 = case_when(
    !is.na(FossilFuelMessenger) ~ FossilFuelMessenger,
    !is.na(CommunityLeaderMessenger) ~ CommunityLeaderMessenger,
    !is.na(EnviroMessenger) ~ EnviroMessenger,
    TRUE ~ NA_character_
  ))

cdr_messenger$support <- factor(cdr_messenger$support, levels = c("Strongly oppose", "Somewhat oppose", "Somewhat support", "Strongly support"))

mess_mod <- polr(support ~ Messenger_party_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr_messenger, Hess=TRUE)
tab_model(mess_mod)

#Check proportional odds assumption
brant(mess_mod)

cdr_messenger$support2 <- factor(cdr_messenger$support2, levels = c("Strongly decrease", "Somewhat decrease", "No difference in my support","Somewhat increase", "Strongly increase"))

mess_mod2 <- polr(support2 ~ Messenger_party_cat+Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr_messenger, Hess=TRUE)
tab_model(mess_mod2)
brant(mess_mod2)

cdr_messenger$MessengerFollowUp_fac <- factor(cdr_messenger$MessengerFollowUp, levels = c("Strongly oppose", "Somewhat oppose","Somewhat support", "Strongly support"))

mess_mod3 <- polr(MessengerFollowUp_fac ~ Messenger_party_cat+Messenger_sector_cat+party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr_messenger, Hess=TRUE)
tab_model(mess_mod3)
brant(mess_mod3)

####Project parameters####

#Public/private ownership preference

cdr_own <- cdr %>% filter(ProjectOwnership != "Haven’t heard enough to say")
cdr_own$publicly_owned_pref <- 0
cdr_own$publicly_owned_pref[cdr_own$ProjectOwnership == 'Publicly owned and operated'] <- 1

table(cdr$ProjectOwnership)/length(cdr$ProjectOwnership)

# Fit exact logistic regression model
proj_own_elr <- logistf(publicly_owned_pref ~ party3*ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data = cdr_own, pl = TRUE)

summary(proj_own_elr)
tab_model(proj_own_elr)

# Calculate variance inflation factors (VIF)
vif_proj_own_elr <- vif(glm(formula = publicly_owned_pref ~ party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum,
                            data = cdr_own,
                            family = binomial))

# Print the VIF results
print(vif_proj_own_elr)


#Test to see who reports not having heard enough to say

cdr$own_dk <- 0
cdr$own_dk[cdr$ProjectOwnership == "Haven’t heard enough to say"] <- 1

# Fit exact logistic regression model
proj_own_elr_dk <- logistf(own_dk ~ party3*ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, pl = TRUE)

summary(proj_own_elr_dk)
tab_model(proj_own_elr_dk)

# Calculate variance inflation factors (VIF)
vif_proj_own_elr_dk <- vif(glm(formula = own_dk ~ Messenger_sector_cat + Messenger_party_cat + party3 + ClimateChange_num + Gender + Age + HouseholdIncomeBracket_num + Race + LocType + education3 + ff_cap_sum,
                            data = cdr,
                            family = binomial))

# Print the VIF results
print(vif_proj_own_elr_dk)



cdr %>% group_by(EconBenefits) %>% summarise(Percentage=n()/nrow(.))
cdr %>% group_by(FossilFuelRole_Nationalization) %>% summarise(Percentage=n()/nrow(.))

cdr$EconBenefits_fac <- factor(cdr$EconBenefits, levels = c("Greatly worsen", "Slightly worsen","No effect",
                                                              "Slightly improve","Greatly improve"))

econ_mod <- polr(EconBenefits_fac ~ Messenger_party_cat+Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(econ_mod)
brant(econ_mod)


#Government funding and ownership

cdr_gov_long <- cdr %>%
  {if ("GovernmentFundingOwnership" %in% names(.)) . else stop("Column GovernmentFundingOwnership not found")} %>%
  mutate(splitColumn = str_split(GovernmentFundingOwnership, ",(?=T)", simplify = FALSE)) %>%
  unnest(splitColumn) %>%
  mutate(splitColumn = if_else(splitColumn == "", NA_character_, splitColumn)) %>%
  dplyr::select(-GovernmentFundingOwnership) %>%
  rename(GovernmentFundingOwnership = splitColumn)

cdr_gov_long <- relevel_columns(cdr_gov_long, GovernmentFundingOwnership = "The government should not provide any funding for CDR")


cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide all the funding for CDR, and projects should be owned and operated by the government"] <- "All"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide all the funding for CDR, but projects should be owned and operated by carbon dioxide removal companies"] <- "All"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide all the funding for CDR, but projects should be owned and operated by local communities"] <- "All"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership =="The government should provide all the funding for CDR, but projects should be owned and operated by fossil fuel companies"] <- "All"

cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide some of the funding for CDR, but projects should be owned and operated by carbon dioxide removal companies"] <- "Some"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership =="The government should provide some of the funding for CDR, but projects should be owned and operated by local communities" ] <- "Some"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide some of the funding for CDR, but projects should be owned and operated by fossil fuel companies"] <- "Some"
cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership == "The government should provide some of the funding for CDR, and projects should be owned and operated by the government"] <- "Some"

cdr_gov_long$GovtFundAmt[cdr_gov_long$GovernmentFundingOwnership =="The government should not provide any funding for CDR"] <- "None"

cdr_gov_long <- relevel_columns(cdr_gov_long, GovtFundAmt = "None")

cdr_gov_long$GovtFundAmt_ord <- factor(cdr_gov_long$GovtFundAmt, levels=c("None", "Some", "All"))

cdr_gov_long <- relevel_columns(cdr_gov_long, GovtFundAmt = "None")


gov_fund_amt_ord <- polr(GovtFundAmt_ord ~ party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr_gov_long, Hess=TRUE)
tab_model(gov_fund_amt_ord)

#Government funding mechanisms
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should not fund carbon dioxide removal"] <- "Don't fund"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should divert funds away from other forms of climate action, including renewable energy production and climate resiliency"] <- "Other priorities"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should tax polluting industries according to the total amount of carbon dioxide they have ever emitted"] <- "Industry"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should tax polluting industries according to how much carbon dioxide they are currently emitting"] <- "Industry"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should tax polluting industries a certain percentage of the total amount of carbon dioxide they have ever emitted"] <- "Industry"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should tax only the wealthiest members of the public"] <- "Public"
cdr$GovtFund3[
  cdr$GovernmentFundingMechanism == "The government should tax the public "] <- "Public"

cdr <- relevel_columns(cdr, GovtFund3 = "Don't fund")

gov_fund_mod <- multinom(GovtFund3~party3*ClimateChange_num + Gender + Age + Income_num + Race + LocType + education3 + ff_cap_sum, data=cdr)
tab_model(gov_fund_mod)
summary(gov_fund_mod)



####Moral hazard####


#Effect of passing CDR law on overall fossil fuels use

cdr$FossilFuelUse_fac <- factor(cdr$FossilFuelUse, levels = c("Greatly increase", "Somewhat increase","No effect",
                                                          "Somewhat decrease","Greatly decrease"))

#model intercepts not significant
ff_mh <- polr(FossilFuelUse_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(ff_mh)

brant(ff_mh)

#Effect of passing CDR law on overall renewable energy use

cdr$RenewablesUse_fac <- factor(cdr$RenewablesUse, levels = c("Greatly decrease", "Somewhat decrease","No effect",
                                                              "Somewhat increase","Greatly increase"))

#model significant!!!
renw_mh <- polr(RenewablesUse_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
#summary(renw_mh)
tab_model(renw_mh)

brant(renw_mh)


#Effect of passing CDR law on overall carbon pollution

cdr$CarbonPollution3[cdr$CarbonPollution=="Greatly increase" | cdr$CarbonPollution=="Somewhat increase"] <- "Increase"
cdr$CarbonPollution3[cdr$CarbonPollution=="Greatly decrease" | cdr$CarbonPollution=="Somewhat decrease"] <- "Decrease"
cdr$CarbonPollution3[cdr$CarbonPollution=="No effect"] <- "No effect"

cdr$CarbonPollution3_fac <- factor(cdr$CarbonPollution3, levels=c("Decrease", "No effect", "Increase"))

cdr$CarbonPollution_fac <- factor(cdr$CarbonPollution, levels = c("Greatly increase", "Somewhat increase","No effect",
                                                              "Somewhat decrease","Greatly decrease"))
#model not significant
carbon_mh <- polr(CarbonPollution_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(carbon_mh)
brant(carbon_mh)


#Effect of passing CDR law on willingness to call representative to urge they decrease fossil fuel use

cdr$FossilFuelUse_CallRep_fac <- factor(cdr$FossilFuelUse_CallRep, levels = c("Not at all likely", "Only a little likely",
                                                                  "Somewhat likely","Very likely"))

#model not significant
ff_callrep <- polr(FossilFuelUse_CallRep_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(ff_callrep)
brant(ff_callrep)


#Effect of passing CDR law on willingness to call representative to urge they increase renewables use

cdr$RenewablesUse_CallRep_fac <- factor(cdr$RenewablesUse_CallRep, levels = c("Not at all likely", "Only a little likely",
                                                                              "Somewhat likely","Very likely"))

#model not significant
renw_callrep <- polr(RenewablesUse_CallRep_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(renw_callrep)
brant(renw_callrep)

#Effect of passing CDR law on willingness to call representative to urge they decrease carbon pollution

cdr$CarbonPollution_CallRep_fac <- factor(cdr$CarbonPollution_CallRep, levels = c("Not at all likely", "Only a little likely",
                                                                              "Somewhat likely","Very likely"))

#not significant
carbon_callrep <- polr(CarbonPollution_CallRep_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(carbon_callrep)
brant(carbon_callrep)


#Effect of passing CDR law on willingness to act to decrease fossil fuel use in one's community

cdr$FossilFuelUse_Community_fac <- factor(cdr$FossilFuelUse_Community, levels = c("Not likely at all", "Only a little likely",
                                                                              "Somewhat likely","Very likely"))
#not significant
ff_comm <- polr(FossilFuelUse_Community_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(ff_comm)
brant(ff_comm)


#Effect of passing CDR law on willingness to act to increase renewables use in one's community

cdr$RenewablesUse_Community_fac <- factor(cdr$RenewablesUse_Community, levels = c("Not likely at all", "Only a little likely",
                                                                              "Somewhat likely","Very likely"))
#not significant
renw_comm <- polr(RenewablesUse_Community_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(renw_comm)
brant(renw_comm)

#Effect of passing CDR law on willingness to act to lower carbon pollution use in one's community

cdr$CarbonPollution_Community_fac <- factor(cdr$CarbonPollution_Community, levels = c("Not likely at all", "Only a little likely","Somewhat likely","Very likely"))

#not significant
carbon_comm <- polr(CarbonPollution_Community_fac ~ Messenger_party_cat+ Messenger_sector_cat*party3+ClimateChange_num+ Gender+Age +Income_num +Race+LocType+education3+ff_cap_sum, data = cdr, Hess=TRUE)
tab_model(carbon_comm)
brant(carbon_comm)


####Role of fossil fuel industry####

cdr_gg2 <- cdr %>% dplyr::select(FossilFuelRole_Nationalization, FossilFuelRole_Untrustworthy, FossilFuelRole_Experience, FossilFuelRole_ProvideEnergy, FossilFuelUse, EconBenefits,Messenger_party_cat,Messenger_sector_cat,party3,ClimateChange_num, Gender,Age,Income_num,Race,LocType,education3,ff_cap_sum) %>% pivot_longer(cols = c(FossilFuelRole_Nationalization, FossilFuelRole_Untrustworthy, FossilFuelRole_Experience, FossilFuelRole_ProvideEnergy),names_to = "category", values_to = "responses")

cdr_gg2$responses <- factor(cdr_gg2$responses,
                           levels = c("Strongly disapprove", "Somewhat disapprove",
                                      "Somewhat approve", "Strongly approve"),
                                    ordered = TRUE)

cdr_gg2 <- cdr_gg2 %>% filter(!is.na(responses))
cdr_gg2 <- relevel_columns(cdr_gg2, category="FossilFuelRole_ProvideEnergy")

cdr_gg_experience <- cdr_gg2 %>% filter(category=='FossilFuelRole_Experience')
cdr_gg_untrustworthy <- cdr_gg2 %>% filter(category=='FossilFuelRole_Untrustworthy')
cdr_gg_energy <- cdr_gg2 %>% filter(category=='FossilFuelRole_ProvideEnergy')
cdr_gg_nation <- cdr_gg2 %>% filter(category=='FossilFuelRole_Nationalization')


mod_exp <- polr(responses ~ party3+ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data = cdr_gg_experience, Hess = TRUE)
tab_model(mod_exp)
brant(mod_exp)

table(cdr_gg_experience$responses)/length(cdr_gg_experience$responses)

mod_untrust <- polr(responses ~ party3+ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data = cdr_gg_untrustworthy, Hess = TRUE)
tab_model(mod_untrust)
brant(mod_untrust)

table(cdr_gg_untrustworthy$responses)/length(cdr_gg_untrustworthy$responses)

mod_energy <- polr(responses ~party3+ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data = cdr_gg_energy, Hess = TRUE)
tab_model(mod_energy)
brant(mod_energy)

table(cdr_gg_energy$responses)/length(cdr_gg_energy$responses)

#Nationalization model significant!!!
mod_nation <- polr(responses ~ party3+ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data = cdr_gg_nation, Hess = TRUE)
tab_model(mod_nation)
brant(mod_nation)

table(cdr_gg_nation$responses)/length(cdr_gg_nation$responses)


response_to_numeric <- function(x) {
  case_when(
    x == "Strongly approve" ~ 4,
    x == "Somewhat approve" ~ 3,
    x == "Somewhat disapprove" ~ 2,
    x == "Strongly disapprove" ~ 1,
    x == "Haven't heard enough to say" ~ NA_real_
  )
}

cdr_dup <- cdr %>% dplyr::select(ResponseId, FossilFuelRole_Untrustworthy, FossilFuelRole_Experience, FossilFuelRole_ProvideEnergy, FossilFuelRole_Nationalization)

cdr_dup <- cdr_dup %>%
  mutate(across(starts_with("FossilFuelRole_"), response_to_numeric)) %>%
  rowwise() %>%
  mutate(FossilFuelRole_favorite = {
    responses <- c(Untrustworthy = FossilFuelRole_Untrustworthy,
                   Experience = FossilFuelRole_Experience,
                   ProvideEnergy = FossilFuelRole_ProvideEnergy,
                   Nationalization = FossilFuelRole_Nationalization)

    if (all(is.na(responses))) {
      "Haven't heard enough to say"
    } else {
      max_response <- max(responses, na.rm = TRUE)
      favorite_roles <- names(responses)[responses == max_response]

      if (length(unique(responses[!is.na(responses)])) == 1 && length(responses[!is.na(responses)]) > 1) {
        "All same"
      } else {
        paste(sort(favorite_roles), collapse = ", ")
      }
    }
  }) %>%
  ungroup()

cdr_dup <- cdr_dup %>% dplyr::select(-FossilFuelRole_Untrustworthy,-FossilFuelRole_Experience,-FossilFuelRole_ProvideEnergy,-FossilFuelRole_Nationalization)

cdr <- left_join(cdr, cdr_dup, by = 'ResponseId')

cdr <- relevel_columns(cdr, FossilFuelRole_favorite = "Haven't heard enough to say")

ff_multi <- multinom(FossilFuelRole_favorite~party3*ClimateChange_num+ Gender+Age +Income_num +Race+education3+ff_cap_sum+LocType, data=cdr)
tab_model(ff_multi)
summary(ff_multi)

table(cdr$FossilFuelRole_favorite)

cdr_sel <- cdr %>% dplyr::select(FossilFuelRole_Untrustworthy,FossilFuelRole_Experience,FossilFuelRole_ProvideEnergy,FossilFuelRole_Nationalization, FossilFuelRole_favorite)
